MODULE pfmm_J_builder

   USE fmm_global_consts
   USE pfmm_global_types
   USE fmm_global_types
   IMPLICIT NONE
   PRIVATE
   ! Public procedures
   PUBLIC :: pfmm_get_J_from_V

CONTAINS

!-------------------------------------------------------------------------------

   SUBROUTINE pfmm_get_J_from_V(scheme,Vff,Tvec_map,LHS,J_matrix,energy)

      IMPLICIT NONE
      TYPE(fmm_scheme),   INTENT(IN)    :: scheme
      REAL(REALK),        INTENT(IN)    :: Vff(:,:)
      TYPE(T_vector_map), INTENT(IN)    :: Tvec_map(:)
      TYPE(raw_mm_data),  INTENT(INOUT) :: LHS
      REAL(REALK),        INTENT(OUT)   :: J_matrix(:,:,:)
      REAL(REALK),        INTENT(OUT)   :: energy

      REAL(REALK)   :: g
      INTEGER(INTK) :: u,v, i,j,k, lm_max

      IF (SIZE(LHS%paras) /= SIZE(Vff,2)) CALL fmm_quit('get_J_from_V')
      IF (SIZE(LHS%qlm_T,2) /= SIZE(Vff,2)) CALL fmm_quit('get_J_from_V:2')
      IF (SIZE(LHS%qlm_T,1) /= SIZE(Vff,1)) CALL fmm_quit('get_J_from_V:3')

      J_matrix = zero
      energy = zero

      ! Build J(mu,nu)^T
      lm_max = MIN(SIZE(LHS%qlm_T,1),SIZE(Vff,1))
      DO u = 1, SIZE(LHS%paras)
         v = LHS%paras(u)%id
         g = DOT_PRODUCT(LHS%qlm_T(:lm_max,v),Vff(:lm_max,v))
         i = LHS%J_indices(v)%i_indx
         j = LHS%J_indices(v)%j_indx
         k = Tvec_map(v)%Tindex
         IF (i > UBOUND(J_matrix,1)) CALL fmm_quit('J_matrix size 1')
         IF (j > UBOUND(J_matrix,2)) CALL fmm_quit('J_matrix size 2')
         IF (k > UBOUND(J_matrix,3)) CALL fmm_quit('J_matrix size 3')
         J_matrix(i,j,k) = J_matrix(i,j,k) + g
      END DO

      ! Get energy after factoring in density to LHS
      CALL pfmm_get_E_from_Vff(scheme,LHS,Vff,Tvec_map,energy)

   END SUBROUTINE pfmm_get_J_from_V

!-------------------------------------------------------------------------------

   SUBROUTINE pfmm_get_E_from_Vff(scheme,LHS,Vff,Tvec_map,energy)

      IMPLICIT NONE
      TYPE(fmm_scheme),   INTENT(IN)  :: scheme
      TYPE(raw_mm_data),  INTENT(IN)  :: LHS
      REAL(REALK),        INTENT(IN)  :: Vff(:,:) 
      TYPE(T_vector_map), INTENT(IN)  :: Tvec_map(:)
      REAL(REALK),        INTENT(OUT) :: energy

      REAL(REALK)   :: g 
      INTEGER(INTK) :: u,v, lm_max

      IF (SIZE(LHS%qlm_T,1) /= SIZE(Vff,1)) CALL fmm_quit('get_E_from_Vff:1')
      IF (SIZE(LHS%qlm_T,2) /= SIZE(Vff,2)) CALL fmm_quit('get_E_from_Vff:2')
      IF (SIZE(LHS%dens) /= SIZE(Vff,2)) CALL fmm_quit('get_E_from_Vff:3')

      lm_max = MIN(SIZE(LHS%qlm_T,1),SIZE(Vff,1))
      energy = zero
      DO u = 1, SIZE(LHS%paras)
         v = LHS%paras(u)%id
         g = DOT_PRODUCT(LHS%qlm_T(:lm_max,v),Vff(:lm_max,v))
         energy = energy + g*LHS%dens(v)
      END DO

   END SUBROUTINE pfmm_get_E_from_Vff

!-------------------------------------------------------------------------------

END MODULE pfmm_J_builder

